library(osmdata)
library(opentripplanner)
library(tidyverse)
library(sf)
library(ggthemes)
library(ggspatial)For this assignment, I am interested in looking at the walksheds and drivesheds to Emergency Shelters in Cambridge. It is important for Cambridge residents to be able to quickly access these shelters during emergency situations.
Emergency Shelters data from: https://www.cambridgema.gov/GIS/gisdatadictionary/Public_Safety/PUBLICSAFETY_EmergencyShelters
Reading the shapefile that was saved to my computer:
emergency_shelters <- st_read(
"data_downloads/PUBLICSAFETY_EmergencyShelters.shp/PUBLICSAFETY_EmergencyShelters.shp")## Reading layer `PUBLICSAFETY_EmergencyShelters' from data source `/Users/jocelyntsai/Desktop/GitHub/ctsai1-vis/data_downloads/PUBLICSAFETY_EmergencyShelters.shp/PUBLICSAFETY_EmergencyShelters.shp' using driver `ESRI Shapefile'
## Simple feature collection with 15 features and 12 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 749626.1 ymin: 2954843 xmax: 767879 ymax: 2968207
## projected CRS: NAD83 / Massachusetts Mainland (ftUS)
The data that I downloaded was originally in NAD 1983 StatePlane Massachusetts Mainland coordinate system, I have to transform it to latitude/longitude first using WGS84 (EPSG:4326) bounds or else I will run into en error when trying to make isochrones later.
The opq function creates a call to the open street map server to download street network data. bbox (bounding box), is like the search bar on google maps (typing in general area name is fine, it doesn’t have to be exact). osmdata_xml will create a file with the information about the street network, and it will save it to the default directory created.
opq(bbox = 'Cambridge MA USA') %>%
add_osm_feature(key = 'highway') %>%
osmdata_xml(file = 'OTP/graphs/default/cambridge_streets.osm')I am transforming my spatial data of street lines onto the NAD83(HARN) / Massachusetts Mainland projected coordinate system from https://spatialreference.org/ref/epsg/2805/ . The add_osm_feature when ‘highway’ is specified takes all the roads in Open Street Map. I then used osmdata_sf to also download data in the sf simple features format so I can plot it.
MA_mainland <- "+proj=lcc +lat_1=42.68333333333333 +lat_2=41.71666666666667 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"
cambridge_street_features <- opq(bbox = 'Cambridge MA USA') %>%
add_osm_feature(key = 'highway') %>%
osmdata_sf()
cambridge_streets <- cambridge_street_features$osm_lines %>%
st_transform(crs = MA_mainland)I then plotted the streets in Cambridge to see what’s included:
Download a Java utility called otp.jar and save it to the OPT folder (only need to run this code chunk once):
## The OTP will be saved to OTP/otp.jar
Build a graph representing street and transit networks:
path_data <- file.path(getwd(), "OTP")
path_otp <- paste(path_data, "otp.jar",sep = "/")
otp_build_graph(otp = path_otp, dir = path_data, memory = 1024) Setting up Open Trip Planner (OTP application opens in the web browser):
Connect to Open Trip Planner:
## Router http://localhost:8080/otp/routers/default exists
I created isochrones for areas within an 8 minute walk and an 8 minute drive of emergency shelters in Cambridge using the otp_isochrone function. I then used the rbind function to combine these sets of polygons (drive, walk) into one dataframe. otp_stop() function stops the Open Trip Planner application from running.
iso_8min_walk <-
otp_isochrone(otpcon = otpcon, fromPlace = emergency_shelters,
mode = "WALK", cutoffSec = 480) %>%
st_transform(crs = MA_mainland) %>%
mutate(mode = "walk")
iso_8min_drive <-
otp_isochrone(otpcon = otpcon, fromPlace = emergency_shelters,
mode = "CAR", cutoffSec = 480) %>%
st_transform(crs = MA_mainland) %>%
mutate(mode = "drive")## Warning in otp_isochrone(otpcon = otpcon, fromPlace = emergency_shelters, :
## Failed to get isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 480},"id":"fid--6db43113_17501b2ddf3_-7ff1"}]}Failed to get
## isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 480},"id":"fid--6db43113_17501b2ddf3_-7fee"}]}Failed to get
## isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 480},"id":"fid--6db43113_17501b2ddf3_-7fe9"}]}Failed to get
## isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 480},"id":"fid--6db43113_17501b2ddf3_-7fe8"}]}
left_side <- st_bbox(iso_all_modes)$xmin
right_side <- st_bbox(iso_all_modes)$xmax
bottom_side <- st_bbox(iso_all_modes)$ymin
top_side <- st_bbox(iso_all_modes)$ymax
ggplot(iso_all_modes) +
geom_sf(data = cambridge_streets, color = "gray") +
geom_sf(aes(fill = mode), alpha = 0.3) +
geom_sf(data = emergency_shelters) +
coord_sf(xlim = c(left_side, right_side),
ylim = c(bottom_side, top_side), expand = FALSE) +
scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
labels = c("By Driving", "By Walking"),
option = "plasma") +
annotation_scale(location = "br")+
annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
theme_map()Typing rosm::osm.types() into the console shows all other basemap options.
ggplot(iso_all_modes) +
annotation_map_tile(zoomin = 0, progress = "none") +
geom_sf(aes(fill = mode), alpha = 0.3) +
geom_sf(data = emergency_shelters) +
coord_sf(xlim = c(left_side, right_side),
ylim = c(bottom_side, top_side), expand = FALSE) +
scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
labels = c("By Driving", "By Walking"),
option = "magma") +
annotation_scale(location = "br")+
annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
theme_map() +
labs(caption = "Basemap Copyright OpenStreetMap contributors")## Loading required namespace: raster
ggplot(iso_all_modes) +
annotation_map_tile(zoomin = 0, type= "hillshade", progress = "none") +
geom_sf(aes(fill = mode), alpha = 0.3) +
geom_sf(data = emergency_shelters) +
coord_sf(xlim = c(left_side, right_side),
ylim = c(bottom_side, top_side), expand = FALSE) +
scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
labels = c("By Driving", "By Walking"),
option = "plasma") +
annotation_scale(location = "br")+
annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
theme_map() +
labs(caption = "Basemap Copyright OpenStreetMap contributors")ggplot(iso_all_modes) +
annotation_map_tile(zoomin = 0, type= "cartodark", progress = "none") +
geom_sf(aes(fill = mode), alpha = 0.3) +
geom_sf(data = emergency_shelters) +
coord_sf(xlim = c(left_side, right_side),
ylim = c(bottom_side, top_side), expand = FALSE) +
scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
labels = c("By Driving", "By Walking")) +
annotation_scale(location = "br")+
annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
theme_map() +
labs(caption = "Basemap Copyright OpenStreetMap contributors")ggplot(iso_all_modes) +
annotation_map_tile(zoomin = 0, type= "cartolight", progress = "none") +
geom_sf(aes(fill = mode), alpha = 0.3) +
geom_sf(data = emergency_shelters) +
coord_sf(xlim = c(left_side, right_side),
ylim = c(bottom_side, top_side), expand = FALSE) +
scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
labels = c("By Driving", "By Walking"),
option = "plasma") +
annotation_scale(location = "br")+
annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
theme_map() +
labs(caption = "Basemap Copyright OpenStreetMap contributors")ggplot(iso_all_modes) +
annotation_map_tile(zoomin = 0, type= "thunderforestoutdoors", progress = "none") +
geom_sf(aes(fill = mode), alpha = 0.3) +
geom_sf(data = emergency_shelters) +
coord_sf(xlim = c(left_side, right_side),
ylim = c(bottom_side, top_side), expand = FALSE) +
scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
labels = c("By Driving", "By Walking"),
option = "plasma") +
annotation_scale(location = "br")+
annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
theme_map() +
labs(caption = "Basemap Copyright OpenStreetMap contributors")ggplot(iso_all_modes) +
annotation_map_tile(zoomin = 0, type= "osmgrayscale", progress = "none") +
geom_sf(aes(fill = mode), alpha = 0.3) +
geom_sf(data = emergency_shelters) +
coord_sf(xlim = c(left_side, right_side),
ylim = c(bottom_side, top_side), expand = FALSE) +
scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
labels = c("By Driving", "By Walking"),
option = "plasma") +
annotation_scale(location = "br")+
annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
theme_map() +
labs(caption = "Basemap Copyright OpenStreetMap contributors")ggplot(iso_all_modes) +
annotation_map_tile(zoomin = 0, type= "thunderforestlandscape", progress = "none") +
geom_sf(aes(fill = mode), alpha = 0.3) +
geom_sf(data = emergency_shelters) +
coord_sf(xlim = c(left_side, right_side),
ylim = c(bottom_side, top_side), expand = FALSE) +
scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
labels = c("By Driving", "By Walking")) +
annotation_scale(location = "br")+
annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
theme_map() +
labs(caption = "Basemap Copyright OpenStreetMap contributors")To calculate the area of each isochrone, I used the st_area() function. The pivot_wider() function then creates a seperate column for each mode of transportation (drive or walk), and each row represents a location with both the drive and walk isochrones instead of having each row in the dataframe just containing either one of the isochrones. For labels, I used breaks/1000000 to convert meters square to kilometer square.
iso_areas <- iso_all_modes %>%
mutate(area = st_area(iso_all_modes)) %>%
st_set_geometry(NULL) %>%
pivot_wider(names_from = mode, values_from = area)
ggplot(iso_areas,
aes(x = as.numeric(walk), y = as.numeric(drive))) +
geom_point() +
scale_x_continuous(name =
"Area within an 8 minute walking distance\nof an emergency shelter (square km)",
breaks = breaks <- seq(10000, 500000, by = 50000),
labels = breaks / 1000000) +
scale_y_continuous(name =
"Area within an 8 minute driving distance\nof an emergency shelter (square km)",
breaks = breaks <- seq(0, 3000000, by = 500000),
labels = breaks / 1000000)+
theme_economist()## Warning: Removed 4 rows containing missing values (geom_point).
ggplot(iso_areas,
aes(x = as.numeric(walk),
y = as.numeric(drive))) +
geom_point(alpha = 0.8, size = 1) +
stat_smooth(color = "lightblue", linetype = 1, size = 0.8) +
scale_x_continuous(name = "Area within an 8 minute walking distance\nof an emergency shelter (square km)",
breaks = breaks <- seq(10000, 500000, by = 50000),
labels = breaks / 1000000) +
scale_y_continuous(name = "Area within an 8 minute driving distance\nof an emergency shelter (square km)",
breaks = breaks <- seq(0, 3000000, by = 500000),
labels = breaks / 1000000)+
coord_polar(theta = "x") +
theme_minimal() ## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing missing values (geom_point).